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Abstract 

We make a detailed numerical study of the spectrum of two Schrodinger operators L± arising 
in the linearization of the supercritical nonlinear Schrodinger equation (NLS) about the standing 
wave, in three dimensions. This study was motivated by a recent result of the second author on 
conditional asymptotic stability of solitary waves in the case of a cubic nonlinearity. Underlying 
the validity of this result is a spectral condition on the operators L±, namely that they have 
no eigenvalues nor resonances in the gap (a region of the positive real axis between zero and 
the continuous spectrum,) which we call the gap property. The present numerical study verifies 
this spectral condition, and further shows that the gap property holds for NLS exponents of the 
form 2/3+1, as long as /3* < < 1, where 

0* = 0.913958905 ± le - 8. 

Our strategy consists of rewriting the original eigenvalue problem via the Birman-Schwinger 
method. From a numerical analysis viewpoint, our main contribution is an efficient quadrature 
rule for the kernel l/\x — y\ in R 3 , i.e., provably spectrally accurate. As a result, we are able 
to give similar accuracy estimates for all our eigenvalue computations. We also propose an 
improvement of the Petviashvili's iteration for the computation of standing wave profiles which 
automatically chooses the radial solution. 

All our numerical experiments are reproducible. The Matlab code can be downloaded from 

http : //www. acm. caltech. edu/~demanet/NLS/ 



1 Introduction 

Suppose that ip(t, x) = e tta2 4>(x) with a/0 and x £ M. d is a standing wave solution of the NLS 

idti) + Aip + |V>| 2/ V = 0, (1) 

where < (3 < -A^ if d > 3 and < < oo if d = 1, 2. Here we assume that (f> = (f>(-, a) is a 
ground state, i.e., 

a 2 (t)- A0 = (/> 2/3+1 , 0>O. 

"The second author was partially supported by the NSF grant DMS-0300081 and a Sloan Fellowship. We would 
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It is known that such (f> exist and that they are radial, smooth, and exponentially decaying, see 
Berestycki, Lions [BerLio] and for uniqueness, Kwong [Kwo]. In one dimension d = 1, these ground 
states are explicitly given as 

nx) = — t (2) 

cosh* 3 ((3x) 

when a = 1 (for other values of a ^ rescale), but in higher dimensions no explicit expression is 
known. From now on, we shall assume that d = 3. 

A much studied question is the stability of these standing waves, both in the orbital (or 
Lyapunov) sense and the asymptotic sense. For the former, see for example Grillakis, Shatah, 
Strauss [GriShaStrl], [GriShaStr2] , Weinstein [Weil], [Wei2], Grillakis [Gri], and for the latter, Bus- 
laev, Perelman [BusPerl], Cuccagna [Cue]. Reviews are in Strauss [Str] and Sulem, Sulem [SulSul]. 

In order to study stability, one generally linearizes around the standing wave. This process 
leads to matrix Schrodinger operators of the form 



H = H + V 



A + a 2 




" -Vi 


-v 2 ' 


A - a 2 _ 


+ 


v 2 





on L 2 (R d ) x L 2 

Conjugating TL by the matrix 



Here, Vy = (J3 + l)</> 2/3 and V 2 = f3<p 2fS . 
1 i 
1 - 



leads to the matrix operator 







iL_ 




with 



L_ = - A + a 2 — ^ 

L+ = - A + a 2 - (2(3 + 1)0 2/3 



The continuous spectrum of both L_ and L + equals 



[a 2 , ool 



Since L-<p = and cf> > 0, it follows 



that zero is a simple eigenvalue and the bottom of the spectrum of L_. Moreover, L+dj(j) = for 
1 < j < 3 so that ker(L + ) C {dj(f> : 1 < j < 3}. In fact, for monomial nonlinearities it is known 
that there is equality here, see Weinstein [W612] 1 , and that there is a unique negative bound state 
ofL+. 

It is known that L_ > implies that the spectrum spec(W) satisfies spec(W) C MUilR and that 
all points of the discrete spectrum other than zero are eigenvalues whose geometric and algebraic 
multiplicities coincide. On the other hand, the zero eigenvalue of Ti has geometric multiplicity four 
and algebraic multiplicity eight provided (3 ^ | , whereas for the L 2 -critical case (3 = | the algebraic 
multiplicity increases to ten. For this see [Weil], [BusPerl] or [RodSchSofl], [ErdSch]. 

In order to carry out a meaningful asymptotic stability analysis it is essential to understand 
the discrete spectrum of TL. The root space at zero was completely described by Weinstein [Wei2]. 
Moreover, it is also well-known that 



spec(W) C R iff (3< 



1 In this paper a restriction (3 < 1 is imposed in d 
full range /3 < 2 by means of Weinstein's arguments 



3, but using Kwong's results [Kwo] allows one to obtain the 
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whereas in the range | < (3 < 2 there is a unique pair of simple complex-conjugate eigenvalues ±ry. 
This latter property reflects itself in the nonlinear theory in the following way: Orbital stability 
holds iff (3 < |, see Berestycki, Cazenave [BerCaz], Weinstein [Wei2], Cazenave, Lions [CazLio], 
and [GriShaStrl], [GriShaStr2] . 

In [Sch] the second author investigated conditional asymptotic stability for the unstable case 
(3 = 1. This analysis depended on the fact that zero is the only eigenvalue of TC in the interval 
[—a 2 , a 2 ] and that the edges ±a 2 are not resonances. The fact that ±a 2 are neither eigenvalues 
nor resonances is the same as requiring that the resolvent (7i — z)~ l remains bounded on suitable 
weighted L 2 (M 3 ) spaces for z close to ±a 2 . 

Using some ideas of Perelman [Per2], it is shown in [Sch] that these properties can be deduced 
from the following properties of L+, L_: Neither L + nor L_ have any eigenvalues in the gap 
(0, a 2 ] and L has no resonance at a 2 . 

In one dimension d = 1, the spectral properties of L_ and L + can be determined completely 
since the generalized eigenfunctions of these operators (more precisely, the Jost solutions) can be 
given explicitly in terms of certain hyper geometric functions, see Fliigge [Flu], Problem 39 on 
page 94. This is due to the special form of the ground state (2). 

Unfortunately, it seems impossible to determine similar properties for the case of three di- 
mensions by means of purely analytical methods. We therefore verify this gap property of L± 
numerically via the Birman-Schwinger method. We will refer to the gap property as the fact that 
L± have no eigenvalues in (0,a 2 ] and no resonance at a 2 . Our main result is as follows. 

Claim 1. There exists a number (3* = .913958905 ± le — 8 so that for all (3* < (3 < 1 the gap 
property holds. 

This statement can also be continued beyond (3=1. We only went up to 1 since (3 = 1 alone 
is needed in [Sch]. In the range (3 < (3*, our numerical analysis shows that the operator L + has 
eigenvalues in the gap (0, 1]. This is perhaps surprising, since it shows that the gap property does 
not hold for the entire L 2 (1R 3 ) super-critical range | < (3 < 2. However, it does hold at (3 = 1. In 
particular, the method of proof from [Sch] does not apply to all (3 > | since it relies on the gap 
property. In contrast, in one dimension d = 1, Krieger and the second author showed that this 
method does apply to the entire super-critical range (3 > 2. In fact, there, the gap property does 
hold for all (3 > 1, see [Flu]. 

In the remainder of this paper, we will be concerned with the description of the numerical 
method, and will study its convergence properties. Note that we did not formulate our main 
result as a theorem because the error bound we give on (3* is based on numerical observations of 
convergence of the method. A fully rigorous justification would involve either an a priori estimate 
for this bound, or a provably reliable a posteriori estimate. For this, we would probably need 
to (1) derive quantitative estimates of regularity and decay for 4>(x) and other functions, with 
explicit values of the constants, and (2) give a full treatment of the propagation of round-off errors 
from a countless number of sources. This type of heroic exercise would contribute little to the 
understanding of the accuracy of the numerical method, so we chose to spare the reader. 
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2 Description of the numerical method 



2.1 The Birman-Schwinger method 

A direct numerical computation of the eigenvalues of L± would be problematic near a 2 , the edge of 
the continuous spectrum. For example, it is unclear if a perceived numerical eigenvalue at 0.99 a 2 
belongs to the gap or if it has escaped from the continuous spectrum. Decay of the corresponding 
eigenfunction might help in making a decision, but this criterion is unacceptable over a truncated 
computational domain. Instead, we will reformulate the problem by the Birman-Schwinger method, 
which we now recall. 

Let H = — A — V, where V > is a bounded potential that decays at infinity. In our case, this 
is H = L± — a 2 1. We would like to filter out positive eigenvalues, so assume Hf = —A 2 / where 
A > and / 6 L 2 . Then g = Uf, where U = \/V satisfies 

g = U(-A + X 2 y 1 Ug. 

In other words, g € L 2 is an eigenfunction of 

K(\) = U(-A + \ 2 ) l U 

with eigenvalue one. Note that K(X) is a compact, positive operator. Conversely, if g € L 2 satisfies 
K{\)g = g, then 

f:=U- 1 g = (-A + X 2 y 1 UgeL 2 

and Hf = —A 2 /. Moreover, the eigenvalues of K(X) are strictly increasing as A —> 0. Hence, we 
conclude that 

#{A : ker(F - A 2 ) + {0}} = #{e > I : ker(K(0) - E) ± {0}} , 

counted with multiplicity. 

Finally, in view of the symmetric resolvent identity, viz. 

(H - z)' 1 = (-A - z)- 1 + (-A - z)- l u\l - U(-A - "V(-A - z)" 1 . 

This shows that the Laurent expansion of (H — z)~ l around z = does not involve negative powers 
of z iff I + U(—A — z)~ 1 U is invertible at z = which is the same as requiring that 

ker{/- Ui-A)-^} = {0} 

because of the Fredholm alternative (assuming that V decays sufficiently fast at infinity to insure 
compactness). In other words, if H has no resonance or eigenvalue at the origin, then K(0) will 
not show an eigenvalue E = 1, and conversely. 

Let us count the eigenvalues {Xj}'jL 1 of K(0) = U^—A)^ 1 !/ (which are all non-negative) in 
decreasing order. Then we arrive at the following conclusion: Let iV be a positive integer. Then 
the operator H has exactly negative eigenvalues and neither an eigenvalue nor a 
resonance at zero iff Ai > . . . > Aat > 1 and Aat+i < 1. 

Note that spectrum of the self-adjoint, compact Birman-Schwinger operator K(0) is discrete 
and robust to numerical perturbations near E = 1, which is the desired numerical effect. Since 
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K(0) is compact, eigenvalues cannot accumulate to E = 1 from below, and are in no way related 
to the continuous spectrum of H. 

In view of the preceding, we therefore need to show the following to justify [Sch]: For (3 = 1, 
the second largest eigenvalue of 2 

47r|x — y\ 

is below one, and the fifth largest eigenvalue of 

yOr)^(y) 



K+(x,y) = (2/3 + 1)- 



Att\x — y\ 



is below one. These properties will then imply the gap property, i.e., L± have no eigenvalues in 
(0, 1] and no resonance at 1. 

2.2 The modified Petviashvili's iteration 

The first step of the numerical method is to find the soliton (f)(x), which is the unique positive, 
radial, decaying solution of 

-A0 + = ^ +1 , (3) 

unique up to translation. As mentioned earlier, <f>{x) is in fact exponentially decaying. A naive 
approach would be to solve a descent equation like 

du . . |2 n 

— = An — u + \u\ y u, 
at 

but, as shown in [BerLioPel], this equation is unstable near the fixed manifold of interest. Instead, 
we will solve the modified Petviashvili 's iteration which reads 

3 

4> n+1 = MZ{I - A)- 1 ^! 2 ^) +^4,^ (4) 

j=l 1 

The initial guess <f>Q can for example be taken as a Gaussian. The choice of constants M n ,i? nj ,7 
and S is crucial for convergence of the iteration, and is given by 

Mn = f^aggfjg, (5 ) 

71,3 ;5£[a,-(i^|2^)]Ade' 

2/3 + 1 

7 = ~2T> (7) 
S = -1/2, (8) 



2 Note that due to the scaling x t-> ax and cj>(x, a) = af> <p(ax) we may assume that a = 1. From now on we set 
<j>(x) = (j>(x,l). 
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where dj = and the hat denoting Fourier transformations. This iteration, without the second 
term, was introduced by Petviashvili in 1976, and convergence was proved recently in [PelSte]. 
The addition of the second term is a minor increment whose purpose is to fix a potential source 
of instability due to numerical discretization and to force the iteration to choose the radial soliton 
(centered at the origin). This will be explained and justified in section 3. 

Numerically, (/ — A) -1 is realized in the Fourier domain, and Fourier transformations are 
implemented via the Fast Fourier Transform (FFT) . Discretization issues are addressed in the next 
section. In practice, the iteration is accelerated via Aitken's method applied pointwise, i.e., 



<j)£(x) = c/) n (x) 



Wi( s ) - 4>n(x)Y 



0n+2(x) - 2<j) n+ i(x) + <j) n (x) ' 



The iteration is stopped when the Euler-Lagrange equation (3) is satisfied up to some very 
small tolerance r in L 2 . The resulting approximation of the soliton will be denoted by <j). 



2.3 Truncation and Discretization 

We now take up the task of computing the eigenvalues of the Birman-Schwinger operators as defined 
in Section 2.1. The first step is to truncate the three-dimensional computational domain to a cube 
of sidelength L centered at the origin and to discretize functions f(x) by evaluating them on the 
regular grid 

Xj = (ji>h>h)jf> ( 9 ) 

with ji, j2, jz integers obeying —N/2 < jk < N/2 — 1. Operators are, in turn, discretized as matrices 
acting on 'vectors' of function samples f(xj). Tools of numerical linear algebra can then be invoked 
to compute the eigenvalues of these matrices. 

Typical values of L and N for which discretizing K± is expected to be reasonably accurate 
are L ~ 20 and N ~ 100. In this context, several vectors of N 3 ~ 10 6 function samples can 
comfortably be stored simultaneously in the memory of a 2005-era computer, but we cannot yet 
afford to manipulate matrices containing iV 6 ~ 10 12 elements. This rules out the possibility of using 
popular approaches such as the QR algorithm, which compute eigenvalues by operating directly on 
the matrix entries. 

Instead, we will resort to a modification of the power method, known as the implicitly restarted 
Arnoldi iteration, which is implemented in Matlab's eigs command [Eigs]. This method has the 
advantage of only requiring applications of the operator to diagonalize, i.e., matrix- vector products. 
For well-conditioned problems, such as the one we are addressing, eigs computes the top eigenvalues 
of the finite matrix up to machine precision, i.e., about 15 decimal digits in Matlab. 

The only remaining issue is then to find a good discretization K± of the Birman-Schwinger 
operators K±, and to quantify the accuracy. Multiplication by 4>(x) to some power will be done 
sample- wise on the grid Xj. Inverting minus the Laplacian in a space of decaying functions over R 3 , 
or equivalently convolving with the fundamental solution G(x) = 47 3i , is a bit more complicated. 
Discretizing G(x) by sampling at Xj is quite inaccurate and is problematic if x = belongs to the 
grid Xj. Dividing by |£| 2 in frequency poses similar difficulties. Instead, for reasons which will be 
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explained in Section 3, we use the following discretization, 



G( Xj )= { 1 Jl V 7 i 2 (10) 



where Si(x) = Jq dt. The discrete (circular) convolution of G(xj) with a vector f(xj) is 

computed efficiently as the multiplication G(xj) f(xj) of their respective FFT, followed by an 
inverse FFT. 

In this context, applying the full Birman-Schwinger operator, say K_, to a vector of samples 
f(xj) consists of the obvious sequence of steps: (1) multiply f(xj) by U{xj) = <j)(xj)P, (2) perform 

an FFT, (3) multiply the result by G(xj), (4) perform an inverse FFT, and finally (5) multiply the 
result by U(xj) again. 

Since the complexity of a one-dimensional FFT in Matlab is O(iVlogiV) operations for most 
values of N (not necessarily a power of two), one application of K± will require O (N 3 log N) 
operations. This is a substantial improvement over the naive matrix-vector product which would 
require 0(N e ) operations. 



3 Convergence analysis 

3.1 The modified Petviashvili's iteration 

In this section we discuss convergence of the iteration (4). The first result in this direction, in the 
case 5 = 0, can be found in ref. [PelSte]. To make this discussion self-contained, we recall their 
argument and apply it to our specific problem. 

Theorem 2. (Pelinovsky, Stepanyants.) Let 4>(x) be the unique radial solution of (3) and 
H}(J$?) denote the subset of all radial functions in iJ 1 ^ 3 ). Consider the iteration (4) with M n 
and 7 given by equations (5), but 5 = 0. Then there exists an open neighborhood N of '(f) in H^(R 3 ), 
in which <f> is the unique fixed point and (4) converges to (p. The iteration is strictly stable in the 
sense that, for all <f>o G M , 

||^ + i - 0||i < (1 - C)\\<f> n -cp\\ u < C < 1, (11) 

where || • ||i is the norm in the Sobolev space H 1 ^ 3 ). 

Proof. Put p = 2/3 + 1. Let us first write down the linearized iteration, about the fixed point 4>, for 
the perturbation w n = 4> n — (j). It reads 3 

w n +i = (1 -p)ja n (f> + p 1 + |g| 2 n , (12) 
3 Throughout this paper, we use the following convention for the Fourier transform: 
HO = J e- lx< f(x) dx, f{x) = ^ j e^f{i) <K. 
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where * denotes convolution, a n comes from the linearization of M n , and is given by 

^ fw n <f>P 

(Up to higher-order terms, we have the asymptotic relation 1 + (1 — p)a n ~ M n .) This formula 
suggests that <p plays a special role for the stability of the linearized iteration. Indeed, we claim 
that we can actually expand w n as 

w n = a n (f) + q n , (13) 

where a n is exactly as defined above, and q n is a remainder. In order to see this, let us introduce 
the operator A = (I — A) _1 H, where H = I — A — p^P~ x . It is easy to check that is it bounded 
and self-adjoint with respect to the H 1 inner product, 

{f,g) = (f,(I-A)g). 

The operator A therefore provides a spectral decomposition of L 2 (R), orthogonal with respect to 
the inner product (•,•). It was noticed in [PelSte], (or by a straightforward extension of their 
argument,) that the spectra of A and H obey 

dim(neg(A)) = dim(neg(W)) = 1, 
nuU(A) = null(W) = 3. 

The first four eigenfunctions of A are precisely <fi and di<f>, with eigenvalues 1 — p < and 
respectively. Equation (13) is just the expansion of w n in this orthogonal system. Since the iterates 
4> n are all radial, the components along di4> are zero. The remainder q n belongs to the space Y p 
defined by 

Y p = {u e L 2 (R 3 ) : (u, <jP) = (u, d t (j?) = 0}. 

In the space Y p , the spectrum of A is strictly positive, bounded from below by the fifth eigenvalue 4 
A5 > and from above by 

(u,Au) . c ju^u) 
X M = sup — — = 1 - p inf — — = 1. 

u (U,U) u (U,U) 

For the last equality we have used the fact that (j>{x) > and 4>{x) — > as \x\ — > 00. 
The recurrence equations for a n and q n can be found from equation (12), 

a n +i = (p-l(p- l))an, 
q n +i = {I ~ A)q n . 

The choice we made for 7 ensures that the component along <j) is immediately put to zero (in the 
linearized iteration.) It is also clear that we have ||<2Vi+i||i < (1 — ^5)||9n||i m i? 1 (K 3 ). In the scope 
of the linearized iteration, equation (11) follows with C = A5. 

Call V the nonlinear operator for the Petviashvili iteration in the case 5 = 0, so that (4) is 
written as n +i = V(4> ri ). We follow [PelSte] and apply the contraction mapping theorem in a 

4 If there is no fifth eigenvalue, then A5 equals the edge of the essential spectrum 
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neighborhood of the pixed point (f>, within the closed subspace H^(M. 3 ). We have already computed 
the Frechet derivative of V at 4>, 

V'{<t>)w = (I - A)P Yp w, (14) 

where Py is the projection onto Y p , orthogonal in H l . Let < e < A5. By continuity of V'(u) as a 
function of u, in the operator H l norm, we can assert that there exists a small open neighborhood 
jV of 4> in which 

||P»||i^i < l-A 5 + e. 

By a standard application of the contraction mapping theorem (see [HutPym] p. 126), the fixed 
point is unique in M and we have the estimate 

\\4>n+l -<t>\\i < (I-A5 + e)\\(t> n - (PWl 

This concludes the proof. □ 

Let us now examine how the above argument generalizes to the case 5^0. The purpose of 
the second term in equation (4) is precisely to put to zero the components along the three basis 
functions di4>, should the initial condition not be radial. This is also useful in the context of the 
numerical realization of V, since numerical round-off errors do not correspond in general to radial 
perturbations. To be precise, the linearized iteration (12) becomes 



hP- 1 



w n +i = (1 - p)ja n <f> + 25 h n^ 3 4> + V 1 + | f | 2 " ' 

j=l,2,3 Kl 

with a n as previously and 



^ J w n di<f 

is the component of w n along di4>, in the natural H 1 inner product. So we have 



+ ^2 b n,jdj4> + q ri 



and the recurrence relation for b n i is 



b n +i,i = (1 + 25)6„,i. 

Choosing 5 = —1/2 as advocated previously will put to zero the components along 3$. 

In practice the modification 5 = —1/2 works very well, see section 4. The analysis of the 
linearized iteration does not pose any difficulty. However, we have been unable to extend the 
argument of theorem 2 to the full nonlinear iteration. The operator V is in general not defined on 
H 1 , because the constants R Hi i involve 3/2 derivatives of <\> n in L 2 . A fortiori, V is not Frechet- 
differentiable for functions in H 1 . Note that the contraction argument cannot work, since the 
soliton <fi is not unique in i? 1 (IR 3 ) - it is unique up to a translation. The question of stability of the 
modified Petviashvili's iteration in the non-radial case is possibly related to the problem of proving 
uniqueness of the radial soliton, which, in itself, is not trivial. 
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3.2 Truncation and Discretization 



In this section we prove spectral convergence of the proposed discretization of the Birman-Schwinger 
operator, which is numerical analysis jargon for (almost) exponential convergence with respect to 
the large discretization parameters L and N/L. Recall that we denote by Xj = {ji,j2,js)jf, j £ 2 3 , 
the nodes of a cubic grid with spacing L/N in all three directions, not necessarily bounded in some 
of the arguments that follow. Obviously, we take L < N . Unless otherwise specified, we are 
considering the operator since K + = (2(3 + l)if_. 

We now formulate our main result. Most of the rest of this section is devoted to its justification. 
We will use the notation (x) = \J\ + |x| 2 . 

Theorem 3. Let K be the Birman-Schwinger operator as above, and K its numerical realization, 
extended to functions of continuous x by sampling and interpolation (see below for details). Let 
5 > be the decay rate of U(x) = 4>{x)P , \U(x)\ < C ■ e -5 l x L Assume that Petviashvili's method 
gives an accurate approximation (j) of the soliton in the sense that, for some e > 0, and denoting 
V = ¥, 



\U( Xj ) - U(xj)\ < C-min( 
Then we have, for all s > 3/2 and f G H 



-S\xj\ 



(Xj)' 

, in exact arithmetic, 



(15) 



e + Le SL/4 + 



(16) 



\\(K-K)f\\ L2 <C s - 

for some constant C s > depending on s. 

Discretization is error-free in the context of the Shannon sampling theorem. Let us introduce 

aited square-inte^ 

3 ) = {ue L 2 (M 3 ) : u(0 = 0, |£| > irN/L}. 



B N / L (R ,i ), the space of band-limited square-integrable functions 



B 



N/L{ 

The hat denotes Fourier transformation. We can then define the sampling operator S, and inter- 
polation operator T, as 

S : B N/L - f, f(x) ^ {f(xj)}. 
T:£ 2 ^ B N/L , {f( Xj )} ^Y, h ( x ~ x j)f( x i)- 

Here h(x) is the interpolating kernel defined by h(t;) = (jf) 3 if |^| < ^jf, and zero otherwise. The 
content of Shannon's sampling theory is that S is an isometry from B N / L to £ 2 , and T is in that 
context both the adjoint and a left inverse for S on its range (hence the interpolation property of 
h.) Note that the properly normalized £ 2 norm is 



\\{f(xj)}j\\p 



( L 



N 



With a slight abuse of notations, let us denote by G the operator of discrete convolution by 
G(xj), defined in equation (10). This particular expression for the weights G(xj) is chosen so that 
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the inversion of the Laplacian is exact on decaying band- limited functions. By 'decaying', we mean 
a member of some weighted L 2 space. More generally, let us introduce weighted Sobolev spaces as 

H^(R 3 ) = {u : (x) m u(x) G H S (R 3 )}, 

and equipped with the norm 



U , m = (X) U 



Of course, L m = H^. 



Lemma 4. We have 

(TGSf)(x) = ±- [ jJ—JMdy, 
4vr J R s \x - y\ 

for all f G B N / L (R 3 ) PI L\( 



Proof. The restriction / G L\ is simply chosen so that the integral makes sense. Since / G B N / L , 
we can express it as 

f(x) = ^2h(x-Xj)f(xj). 

3 

Substituting in the integral, we get 

(-A)-V(^) = -L I T^—My-x k )dyf(x k ) 

The quadrature weights w(j, k) can be computed explicitly. By Parseval, 

In the event Xj / i^, we can express this integral in a spherical coordinate frame whose z-axis is 
aligned with Xj — x k . It becomes 

1 / T \ 3 r wN / L r w i 

«,(,", fc) = _^-J j( r 2 dr j o smede- 2 e-^-^™<>, 

dr / e^-^dx, 



1 r 



ivy (2vr)2y y_! 



NJ (2tt) 2 J q \xj-x k \r 
L\ 5 1 1 „. f irN\xj — Xk\ 



iV / 2tt 2 \xj — x k \ \ L 

Si(x) is the special function sini/t dt. Thus we have w(j,k) = G{xj — x k ), as in equation (10). 
The special case Xj = x k can be handled separately in a similar fashion, and one can check that it 
corresponds to the limit Xj — x k —>■ in the general expression. 
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We have just proved that GSf = S(— A" 1 )/. Now the output (— A) -1 / is not in general square- 
integrable, see also Lemma 8 below, but it is still band-limited. The interpolation operator T is 
still the inverse of sampling for band-limited functions with some growth, so we obtain the desired 
conclusion. □ 

Remark. It is interesting to notice that the choice we made for h, the indicator of a ball, is 
not standard. The most natural frequency window associated to a Cartesian grid is the indicator 
of the cube [—ttN/L,ttN/L} 3 . This choice of window for the definition of B^ L , the space of I? 
band- limited functions, would have made the sampling operator S not only an isometry, but a 
unitary map from B^j^ onto £ 2 . In our context, with the spherical window, S is only unitary from 
Bn/l to a subset of £ 2 , namely the range SB^/^. However, we do not believe the cubic window 
would have allowed us to formulate a closed- form expression for the weights G(xj). 

The functions we have to discretize, like the soliton (f>(x), are unfortunately not band-limited. 
The smoother the function, the more accurate its sampling, however. The following elementary 
lemma establishes this property. 

Lemma 5. For all s > 3/2 and f G H S (R 3 ), 

\\f-TSf\\ L ,<C s (^ 



L 

The constant C s is a decreasing function of s £ (3/2, oo). 

Proof. Sampling in x corresponds to periodizing in £, so the Fourier transform of TSf — f is 

f(o = /(£ - 2 ™f )x\&<«n/l(z) - no- 



We get two contributions, call them (I) and (II), in 



1*11! = / I E /(£ - 27Tm 

J\Z\<*N/L m ^ Q 
+ / \f(0\ 2 dt. 

J\£\>ttN/L 



L ' 



di 



l\£\>nN/L 

The second integral (II) is easily bounded by 



\£\>7rN/L 

The first integral (I) can be expressed as 



L 



-2s 



(II) < sup (C)- 2s • ll/ll? < vr- 2 ' i? A ] 



s" 



(^) - X] E Im ' m> 
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where 



I,,,.,,,'- I |/(e-27rm^)/(e-27rm'^)|, 



/ 

JBnin 



B (nN/L) L L 

N N 

< sup (£ — 2itm—) s ■ sup (£ — 2irm'—) 

\$\<irN/L L \£\<itN/L L 



[' ~ N N N N 

x / - 27rm-)|(e - 2vrm-) s • - 2irm'-)\(Z - 2m'-) s d£. 

JB (ttN/L) l l l l 



Obviously sup^| <7rA r / / i (^ — 2irm^) s < C • (|m|^) s , and the integral can be bounded by Cauchy- 
Schwarz. We have thus separated 



Im.m' — C ' JmJm' 



N \ -2s 

T 



where 



J m = \m\-\ [ \f(0\ 2 (0 2s dC 

y J|?+27rmf \<wN/L 

Each point £ in the frequency domain is covered by at most one ball B 27Tm N_(7rN/L). Consequently, 
Cauchy-Schwarz for sequences gives 



Y. J ™<c- /£ h- 2s • Jf R3 \mw(0 2s 



m^O y m^O 

The sum over m converges only when s > 3/2 in three dimensions. This governs the dependence 
of the constant C s on s. Each J m gives one factor ||/|| s , which completes the argument. □ 

Notice that the above result is sharp with respect to the range of s for which it is valid. 
The Sobolev embedding from H s into continuous functions is only valid when s > 3/2 in three 
dimensions, and it does not in general make sense to sample a discontinuous function. 

Let us note in passing that Lemma 5 generalizes without difficulty to weighted norms. 

Lemma 6. For all s > 3/2, m G Z, and f G H° n (M. 3 ), 

\\f-TSf\\ L?n <C s (^ 



Proof. It is sufficient to exhibit a bandlimited multiplier £ m (x) equivalent to (x) m in the sense that 

C 1 (x) m <£ m (x)<C 2 (x) m , 
for then weighted norms can be expressed with £ m instead, and 

U{f -TSf) = {Id-TS)Uf. 
Then the conclusion would follow from an application of Lemma 5 to £ m f. 
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Such a function £ m (x) can be constructed by convolving (x) m by some appropriate nonnegative, 
band-limited kernel. For example, 



x 



u(x)={x) m *n( smc (^)) 2H+4 



a 

3=1 

will do for m / 0, provided a > is chosen so that the band limit of the kernel is compatible with 
the sampling of S. □ 

We will frequently need to switch from discrete to continuous norms, using the following result. 

Corollary 7. Let f G H S (R 3 ), for s > 3/2. Then 

H/O* <C S - \\f\\ H .. 

Proof. Decompose TSf as / + (TSf — /), and estimate in L 2 . By Lemma 5 and the sampling 
theorem, 

<II/IIl* + c, •(£) S \\f\\ Hs <c \\f\\ Hs . 

□ 

Let us collect some more background results. We have already hinted at the fact that (— A) -1 , 
defined as the convolution with G(x), is not bounded on L 2 , or on any H s for that matter. But 
it is bounded between weighted Sobolev spaces. Let us recall the following classical result from 
[ JenKat] . 

Lemma 8. (Jensen-Kato) The kernel G(x — y) = , . _ . maps boundedly H' 1 ^ 3 ) toH^_ m ,(M. 3 ), 
and is in addition Eilbert- Schmidt from to L 2 _ m ,, provided 

m,m' > — and m + m! > 2. 

A similar property holds for the discretized kernel G. Weigthed discrete £ 2 spaces, relative to 
the grid defined as 

e 2 m = { Uj : ( Xj ) m Uj e e 2 }. 

Lemma 9. The kernel G(xj — x^) defined in equation (10) maps boundedly l? m to i 2 _ m i provided 

m,m' > — and m + m! > 2. 

The operator norm of G is uniform in N/ L > 1 . 

Proof. The Hilbert-Schmidt norm of G between and £ 2 _ m , is 

iiGii2r S = E - x k )\ 2 (x 3 r 2m '( Xk r 2m . 
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The diagonal part D of the sum, for j = k, is 

d = — ( L\ A Y(i + \i-?r {m+m,) 

(2tt) 2 \n) 2^^^\J N \ > 
and can be compared to the integral 

D<C-^- [ (1 + |x| 2 )- (m+m,) dx, 
which converges when m + m' > 3/2. The off-diagonal part OD, for j / k, is bounded by 



od<[-) y j 1 T (i + \j^y m '(i + \k-\ 2 )- m 

(because | < 2.) Its continuous counterpart is 

ODKC-! [ , 1 t J l + \x\ 2 y m (l + \y\ 2 r m ' dxdy. 

The double integral is the Hilbert-Schmidt norm of G, squared, between L 2 ^ and L 2 _ m/ , and is 
bounded by lemma 8. □ 

We can now turn to the proof of the main result, theorem 3. The Birman-Schwinger operator 
is defined as 

K = U(x)(-A- 1 )U{x), U{x)=(j){xf. 
The approximation K is defined as 

K = TU{ Xj )G c U( Xj )S, 

where the subscript c in G c indicates that the convolution is not over the whole infinite grid 
L/N x Z 3 , but is a circular convolution over the finite cubic array (-L/2 : L/N : (L/2 — L/N)) 3 
(in Matlab notation.) Let us denote this bounded grid by a N . 

Proof of theorem 3. Let us divide the proof into four successive approximation steps, using the 
triangle inequality. 

1. Let us first show that 

(I) = \\TU(x j )(G c -G)(U(x j )Sf)\\ L 2 

is adequately small in the sense of theorem 3. By Shannon's sampling theorem (ST is the 
orthogonal projector onto RanS 1 ,) we can rewrite 

(I) < \\U( Xj )(G c - G)(U( Xj )f( Xj ))y- 

We denote the operation of folding back a grid point by periodicity onto the grid O n as 
follows: 

[ Xj ] = ( Xj + (1, 1, l)£)mod L - (1, 1, 1)£. 



15 



The discrepancy between the two types of convolution is 

(G c - G)g{xj) = g(x k )G{xj - x k ) - ^ g{x k )G{\xj - x k ]). 

kez 3 x k en N 

Let us introduce the intermediate quantity 

G b g(xj) = 9{xk)G{xj - x k ), 

where the subscript b stands for 'bounded convolution'. The first contribution is 

(I A ) = \\U(x j )(G b -G)(U(x j )f(x j ))\\ P , 

= \\U(xj) G(xj - x k )U(x k )f(x k )\\ e 2, 

x k £u N 

< sup|f/(x i )(x j ) 2 | • ||G||^2 2 • sup \U(xj)(xj) \ ■ \\f(xj)\\fi. 

The first factor is obviously bounded by equation (15), the second factor is bounded by 
Lemma 9, the third factor is less than e + C ■ Le~ SL by equation (15), and the fourth factor 
is less than C • ||/||h* by Corollary 7. 
The second contribution is 

(I B ) = \\U(x j )(G c -G b )(U(x j )f(x j ))\\ P . 

The kernels G c and G b differ only when Xj ^ N: and in that case we have the estimate 



\G{xj) - G{[ Xj ])\ < T Xx^n N {j,k). 



Therefore 



(Ib)< 



U{ Xj ) Y, U(x k )f(x k ) 



x k en N 

x k £xj+U N 



1 

< — 

- L 



U{xj) 



\ 



x k en N 
x k £xj+n N 



■ \\f(xj)y. 



The quantity underneath the square root can be bounded, up to a multiplicative constant, 
by 



L 



x£B Xj (L/2) 



e- 2S ^dx< [ e- 2& Wdx. 

JxfB {L/2-\xj\) 

< C .(^-\ Xj \) 2 e- 2S ^-^. 
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With this bound, (/#) becomes 
C 



(ft) < L 



\ / f en »r 



e- s \*\£-\x\)e- s &-Wdx- \\f{ Xj )\\p, 

B (L) J 



<CLe-' L ^\\f(x j )y. 

As before, the I 2 norm of f(xj) can be bounded, up to a constant, by the H 8 norm of / 
(Corollary 7). This shows that (I a) and (Ib) are both within the bounds of equation (16). 

2. Let us now study the difference 

(II) = \\TU(x j )G(U(x j )Sf - SU(x)f)\\ L2 

We have already observed that TU(xj) is bounded between £^_ 2 an d L 2 • We also know from 
Lemma 9 that G is bounded from £ 2 to £ 2 _ 2 - Hence, it suffices to show that the following 
quantity is adequately small: 

\\(U(x 3 )S-SU(x))f\\ 2 } = (L\'Y^{x 3 ) 2 \U(x 3 )-U(x 3 )\ 2 \f(x 3 )\ 2 . 

3 

By equation (15) and Corollary 7, we can bound this expression by Ce 2 ||/||^ s . 

3. The third contribution 

{III) = \\(TU(xj) - U(x)T)GSU(x)f\\ L , 

can be bounded analogously. We know that SU maps £ 2 to l\ boundedly, and G maps l\ to 
£ 2 _ 1 boundedly. It remains to show that 

\\(TU( Xj )-U(x)T)g\\ L i 

is small in proportion to g G l 2 _ x . By the sampling theorem, this is also 

\\(U( Xj ) - SU(x)T)g( Xj )y = \\(U( Xj ) - Uix^g^y 

Again, equation (15) allows to bound this quantity by Ce||<7(xj)||^2 . 

4. The last, and perhaps most important contribution, is 

(IV) = \\U(TGS + A- 1 )Uf\\ L 2. 

Obviously, multiplication by U(x) is bounded from H s to H 2 , as well as from Ht 2 to H s , for 
all s > 0. So it suffices to show that 

{{(TGS + A-^gll^KCs-f^) 8 -\\g\\H S , 
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for all s > 3/2. It is key to notice that 

TGSg = -A^TSg, 

which follows from applying Lemma 4 to TSg. Since by lemma 8, A -1 is bounded between 
L\ and L 2 _ 2 , it is enough to check that 

\\TSg-g\\ Ll <C s -[^\ ■ \\g\\ H .. 
This is precisely the content of lemma 6. The proof is complete. 

□ 

It is very likely that Theorem 3 could actually be formulated with exponential decay in N/L, 
but we have been unable to extend the argument. In order to do so, it would be necessary to prove 
analyticity of the soliton <p(x), which we believe is still an open problem in three dimensions. Our 
numerical experiments strongly support that conjecture (see Section 4). 



3.3 Accuracy of eigenvalues 

Let us remark that accuracy of the discretization depends crucially on the smoothness of the 
function to which the operator is applied. This follows not only from sampling issues, but also 
from the choice of quadrature weigths G(xj) used to invert the laplacian. As a consequence, 
only the eigenvalues of K corresponding to very smooth eigenfunctions will be computed to high 
accuracy. 

Since K is compact and self-adjoint on L 2 (]R 3 ), let us write its spectral decomposition as 



Its discrete counterpart K is also self- adjoint and compact (of finite rank), but only on the space 
B N / L . We have 

Ke.j = XjSj, (ej,e/c) = 5jk- 

Any / G B N / L can therefore be expanded as / = Ylj{f> 

The following result about accuracy of eigenvalues is mostly a consequence of theorem 3. Let 
£l,n denote the small factor in the right-hand side of equation (16), namely 



£l,n = C s ■ 



< L, -< )L ! + ( ^ 



Corollary 10. If A, is simple, then 



£L,N\\ej\\H s 
\(ej,ej)\ 



(17) 



and, for €l,n sufficiently small, 



e j\\L 2 S 



< 2 e L,N\\ej\\Hs 
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where dj = min^j | Xj — X k | . 

If Xj has multiplicity p > 1, denote by Uj the orthoprojector onto the j-th eigenspace, and Tlj 
its numerical counterpart. Then, for all 1 < m < p, 



~ e L,N\jJ2n=l \\ e j,n\\ 2 Hs 
\Xj - Xj,m\ < |ppj-x |i , (18) 



and, for e^jv sufficiently small, 

e L,N^YZ=l \\ e j,n\\ 2 Hs 



Uj-UjWhs < 2- 



dj 



The above norm is the Hilbert- Schmidt (HS) norm. 

Proof. The proof is loosely related to the argument behind Gershgorin's circle theorem. 

• Take Xj simple. Let II = TS be the orthogonal projection from L? to B N / L . We can then 
expand 

Ilej = ^2 9j tk e k , 9j tk = {ej,e k ). 
k 

Consider the relation 

h '<j ■ '•/• ( 19 ) 

where by theorem 3, the remainder rj = (K — K)ej obeys 

IkjlU 2 < £L,N\\ e j\\H s ■ (20) 
We can project equation (19) onto B N / L and expand it on the basis e k , which gives 

Oj, k Xk = Oj,kXj + (rj,e k ). (21) 

For k = j, we can bound 



which is exactly (17) after using the estimate (20) on rj. 

In order to obtain the estimate for the eigenfunctions, we should instead consider 

Kej = Xjij+fj, 

with 1 1 fj \\ L 2 < €L,N\\ej \\h s , and this time expand it on the basis e k , 

8j,kXk = 0j t kXj + (fj, ejfc). (22) 

For k = j, it would give the same estimate as (17), but with ej substituted for ej. For k / j, 
we can rewrite equation (22) as 

\e hk \ < Ji^I (23) 

\Xj - x k \ 
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We are not quite finished since the eigevalue gap, at the denominator, is measured with \j 
instead of Aj. Put Sj = \Xj — Xj\. Let us also introduce dj = min^j |Aj — A&|, and observe 
that dj > dj — 5j. Then, squaring (23) and summing over k / j, we get 



r .112 



k+j a j 

and, as a result, 

|%,/>l-fep^ (24) 
d 3 

Our estimate for Sj, equation (17), can therefore be improved to 

2 II ||2 
<5j < o ~ ,„ ■ 



1 



It is not hard to see that, as 6l,n gets small, Sj, defined here implicitly, decreases to zero. 
Take L, N and l/e sufficiently large so that Sj < dj/4. Going back to equation (24), we 
obtain 

|%,/>l-2^. 
The estimate for eigenvectors follows since 

\\ej-ej\\ 2 = 2-26jj<2-2e]j. 

Assume now that Aj has multiplicity p > 1. The previous argument goes through, with the 
following modifications. The change of basis coefficients are now 6j n ^ = (ej n ,ek), where 
n = 1, . . . ,p. Equation (21) becomes 

Ojn,k^k = 8jn,k^j + (rj, n ,ek)- (25) 

Let us study the discrepancy between Aj and Xj m , for some m. We need to take a linear 
combination of equation (25) with the weights 



VEn \ 6jn,jm\ 2 

We then obtain the inequality 

Tn a n\( r j,nj &jr 



\Xj - Xj rn \ < 



n,jm\ 



Observe that ||ILjej m ||£2 = y/Yln \ 9j n! j m \ 2 , and the numerator can be bounded by Cauchy- 
Schwarz (in n) to yield (18). 
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As for the eigenfunctions, Hilbert-Schmidt norms and trace inner products are to be substi- 
tuted for their L 2 scalar counterpart. More precisely, 

l|n J -n J ||^5 = Tr(n,-n,) 2 , 

= TrlL, + Trflj - 2TrII j n j , 

= 2p- 2Trn j n j , 

so we need to find a bound on the latter quantity. Observe first that 

v 

™A = E E i%«^i 2 - 

n m=l 

The change of basis coefficients are normalized in the sense that, for all n = 1, . . . ,p, 

1 = |#jfn,fc| = \ @jn,jm\ + |#jfn,/c| • 
k m k^j 

Summing over n, we get 

p = TrnA'+EEl^| 2 - ( 26 ) 

n k^j 

The bound on 9j n ^ for k ^ j can be obtained as above, and is 

Wjn,k\ < 7~ — ■ (27) 

It is natural to define c/j = minfc in |Aj n — A^|. The way to bound dj from below by, say, 3dj/4 
is very analogous to the simple case, and based on some adequate control of the size of the 
denominator in equation (18). This can be obtained from 

1 = 1 1 Sjjj! 1 1 2 + ^ ^ \dk,jm\ j 
k^jn 

and 

I 2 < V \( rk > l 3m)? = \(ek,r jm )\ 2 < \\r jm \\ 2 L 2 

(we have used K = K* and K = K* for band- limited functions.) Eventually, we use theorem 
(3) one more time to bound 

||^'n||z,2 < eL,N\\ejn\\H>>- (28) 

The desired estimate is obtained by combining the intermediate inequalities (26) to (28). This 
concludes the proof. 

□ 

Note that, in the formulation of the corollary or in its justification, it is nowhere necessary to 
obtain bounds on errors in computing other eigenvalues A& ^ Xj. This is the sense in which our 
argument is close to Gershgorin's theorem: we only need one row of the matrix (ej, Ke^), namely, 
the j-th row. 
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4 Numerical results and discussion 



Our numerical approximation to the soliton 4>(x), obtained by Petviashvili's iteration, is plotted in 
log-scale in Figure 1. For this experiment, we have taken (3 = 1, a cube of sidelength L = 20 and a 
grid of size N = 200. Notice the apparent exponential decay both in space and frequency. Figure 
2 checks the convergence of Petviashvili's iteration, namely that M n — > 1 and that — A(j) n -\- (p n — 
<lh? +l —> as n — > oo. The constants R n j are not plotted; in all our tests they are at machine 
accuracy, or below. In other words, the soliton is as radial as allowed by the grid. 



Soliton Fourier transform of the soliton 




Figure 1: Left: the soliton <j>(r) in space. Right: the soliton </>(|£|) in frequency. Both depend only 
on the radial coordinate. 

In the above setting we observe nice convergence of the error, as measured by different indicators, 
to machine accuracy. This is not always the case. For values of L too big in comparison to 
N (typically, for L of the order of N/4 and above,) the Petviashvili iteration possibly reaches 
acceptable error levels, but then diverges away from the fixed point. A careful inspection of the 
remainder — A(f> n + 4> n — \<f>n \4>n indicates that this might be due the fifth eigenvalue A5 of the 
linearized iteration operator A becoming negative, in the notations of Section 3.1. 

Next, we show a plot of the two, resp. five largest eigenvalues of the Birman-Schwinger operator 
resp. K + , as a function of (3 in the range [|, 1]. As mentioned earlier, \2(K) is less than 1 
for all values of (3, but there exists a number /3* below which \$(K + ) > 1. This is the signature 
of at least one eigenvalue in the gap of L + , for 2/3 < (3 < (3* (inspection of A6 reveals that there 
is only one eigenvalue in the gap.) For this experiment, we used L = 15 and N = 60. Note that, 
in both cases, A2 = A3 = A4 is triple, and A5 is simple, so we only show at most 3 curves. Our 
numerical implementation correctly picks up the multiplicity, exactly (to all 16 digits), and in all 
the cases that we have tried. 

An accurate computation of the numerical value of the exceptional exponent /?* requires higher 
values of L and N. On a 2005 standard desktop, we have tried L = 25 and N = 200. We then 
obtain [3* by interpolation of As(K + ) for different closeby values of /?, see table 1. The confidence 
on /?*, which we estimate to be about 8 digits, is directly related to the level of accuracy of \r,(K + ). 
The latter is determined by inspection of convergence as L and N increase. So the bounds given 
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Convergence of Petviashvili's iteration 

10 6 | , , , , , , r- 




5 10 15 20 25 30 35 40 45 



Number of iterations 



Figure 2: Convergence of Petviashvili's iteration. The x-axis is the number of iteration. Dotted 
line: 1 — M n , where M n is the Petviashvili constant. Dash-dotted line: Euler-Lagrange remainder, 
|| — A(f) n + 4>n — 0n' 9+1 ||/||0n||, with norms in L 2 . Solid line: Same remainder, but for the Aitken 
iterate 4>^. 



in the main claim, in the introduction, are not rigorous, but merely serve as an indication of how 
accurate we believe our algorithm is. 



0.91395850 1.00000016477 

0.91395875 1.00000006304 

0.91395900 0.99999996130 

0.91395925 0.99999985957 



Table 1: Fifth eigenvalue of K+, as a function of [3 near Here L = 25 and N = 200. Each value 
took about one day to obtain. Cubic interpolation reveals /3* ~ 0.913958905 ± le — 8. 

At this point the reader might wonder if, instead of a full three-dimensional simulation, there 
exists a computational strategy involving only the radial coordinate to compute both the soliton 
and the eigenvalues of the Birman-Schwinger operators. We believe the answer is positive, but will 
involve significantly different ideas from the ones presented in this paper. In particular, spectral 
accuracy will be more difficult to obtain. We leave the problem of determining more digits of the 
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Figure 3: Left: Two largest eigenvalues of K_, as a function of (3 in the range [|, 1]. Right: Five 
largest eigenvalues of K + , as a function of (5. In both cases, the largest eigenvalue is simple, the 
second to fourth eigenvalues are one triplet, and the fifth eigenvalue is simple. The value one is 
indicated by the dotted line. 



constant /3*, likely through the use of a one-dimensional method, as a challenge to the interested 
reader. 

Finally, the Matlab code we used to generate the figures and compute an estimation of /?* can 
be freely downloaded from 

http : //www. acm. caltech.edu/~demanet/NLS/ . 
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